#########################################################
# Application
# CBQ estimation: 
# Replication of Alvarez and Nagler (1995)
# Q1
#########################################################
# Prepare dataset
  
  library(mlogit)
  source("replication_AN/cbq_function.R")
  dat <- read.table("replication_AN/nes9212r.asc")
  varnames <- c("preschc", # choice: 1 Bush 2 Clinton 3 Perot
      "educ", # Respondents education
      "east", # Region (East)
      "south", # Region (South)
      "west", # Region (West)
      "women", # Gender (Female)
      "respfinp", # Felt personal finance were worse
      "natlec", # Felt national economy was worse
      "bclibdis", # Ideological Distance
      "gblibdis", # Ideological Distance
      "rplibdis", # Ideological Distance
      "dem", # Democrat
      "rep", # Republican
      "respgjob", # Oppose government jobs
      "resphlth", # Oppose government health care
      "respblk", # Oppose government minority assistance
      "respab", # Abortion
      "termlim", # Term limits
      "age1829", # Age: 18-29
      "age3044", # Age: 30-44
      "age4559", # Age: 45-59
      "newvoter", # New or returning voter
      "deficit") # Felt deficit was a major problem
  names(dat) <- varnames
  dat$choice <- "Bush"
  dat$choice[which(dat$preschc==2)] <- "Clinton"
  dat$choice[which(dat$preschc==3)] <- "Perot"
  names(dat)[9] <- "dist.Clinton"
  names(dat)[10] <- "dist.Bush"
  names(dat)[11] <- "dist.Perot"
  vote_dat <- mlogit.data(dat, varying = c(9:11), shape = "wide", choice = "choice")
  vote_dat$bush <- ifelse(vote_dat$alt=="Bush",1,0)
  vote_dat$clinton <- ifelse(vote_dat$alt=="Clinton",1,0)
  vote_dat$perot <- ifelse(vote_dat$alt=="Perot",1,0)
  vote_dat$respfinp.bush <- vote_dat$respfinp * vote_dat$bush
  vote_dat$respfinp.clinton <- vote_dat$respfinp * vote_dat$clinton
  vote_dat$natlec.bush <- vote_dat$natlec * vote_dat$bush
  vote_dat$natlec.clinton <- vote_dat$natlec * vote_dat$clinton
  vote_dat$respgjob.bush <- vote_dat$respgjob * vote_dat$bush
  vote_dat$respgjob.clinton <- vote_dat$respgjob * vote_dat$clinton
  vote_dat$resphlth.bush <- vote_dat$resphlth * vote_dat$bush
  vote_dat$resphlth.clinton <- vote_dat$resphlth * vote_dat$clinton
  vote_dat$respblk.bush <- vote_dat$respblk * vote_dat$bush
  vote_dat$respblk.clinton <- vote_dat$respblk * vote_dat$clinton
  vote_dat$respab.bush <- vote_dat$respab * vote_dat$bush
  vote_dat$respab.clinton <- vote_dat$respab * vote_dat$clinton
  vote_dat$east.bush <- vote_dat$east * vote_dat$bush
  vote_dat$east.clinton <- vote_dat$east * vote_dat$clinton
  vote_dat$south.bush <- vote_dat$south * vote_dat$bush
  vote_dat$south.clinton <- vote_dat$south * vote_dat$clinton
  vote_dat$west.bush <- vote_dat$west * vote_dat$bush
  vote_dat$west.clinton <- vote_dat$west * vote_dat$clinton
  vote_dat$newvoter.bush <- vote_dat$newvoter * vote_dat$bush
  vote_dat$newvoter.clinton <- vote_dat$newvoter * vote_dat$clinton
  vote_dat$termlim.bush <- vote_dat$termlim * vote_dat$bush
  vote_dat$termlim.clinton <- vote_dat$termlim * vote_dat$clinton
  vote_dat$deficit.bush <- vote_dat$deficit * vote_dat$bush
  vote_dat$deficit.clinton <- vote_dat$deficit * vote_dat$clinton
  vote_dat$dem.bush <- vote_dat$dem * vote_dat$bush
  vote_dat$dem.clinton <- vote_dat$dem * vote_dat$clinton
  vote_dat$rep.bush <- vote_dat$rep * vote_dat$bush
  vote_dat$rep.clinton <- vote_dat$rep * vote_dat$clinton
  vote_dat$women.bush <- vote_dat$women * vote_dat$bush
  vote_dat$women.clinton <- vote_dat$women * vote_dat$clinton
  vote_dat$educ.bush <- vote_dat$educ * vote_dat$bush
  vote_dat$educ.clinton <- vote_dat$educ * vote_dat$clinton
  vote_dat$age1829.bush <- vote_dat$age1829 * vote_dat$bush
  vote_dat$age1829.clinton <- vote_dat$age1829 * vote_dat$clinton
  vote_dat$age3044.bush <- vote_dat$age3044 * vote_dat$bush
  vote_dat$age3044.clinton <- vote_dat$age3044 * vote_dat$clinton
  vote_dat$age4559.bush <- vote_dat$age4559 * vote_dat$bush
  vote_dat$age4559.clinton <- vote_dat$age4559 * vote_dat$clinton
  vote_dat$Y <- ifelse(vote_dat$choice=="TRUE", 1, -1)
  vote_dat <- vote_dat[order(vote_dat$chid,vote_dat$Y),]
  
  Y <- vote_dat$Y
  ind <- vote_dat$chid
  N_indx <- length(unique(ind))
  N <- length(Y)
  X <- as.matrix(subset(vote_dat,select = c(
    dist , respfinp.bush , natlec.bush , respgjob.bush , resphlth.bush , 
    respblk.bush , respab.bush , east.bush , south.bush , west.bush , newvoter.bush , termlim.bush ,
    deficit.bush , dem.bush , rep.bush , women.bush , educ.bush , age1829.bush , age3044.bush , age4559.bush
    , respfinp.clinton , natlec.clinton , respgjob.clinton , resphlth.clinton , 
    respblk.clinton , respab.clinton , east.clinton , south.clinton , west.clinton , newvoter.clinton , termlim.clinton ,
    deficit.clinton , dem.clinton , rep.clinton , women.clinton , educ.clinton , age1829.clinton , age3044.clinton , age4559.clinton
  )))
  n_covariate <- dim(X)[2]
  
  cbq_fit <- cbq(N = N,
                n_covariate =  n_covariate,
                X = X,
                Y = Y,
                N_indx = N_indx,
                ind = ind,
                qtl = 1,
                nchain=2,
                niter = 1000
  )
  
  save(cbq_fit,file="replication_AN/cbq_vote_q1.RData")
      
